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Abstract 

We present a method for solving a class of initial valued, coupled, non-linear differential equa- 
tions with 'moving singularities' subject to some subsidiary conditions. We show that this type 
of singularities can be adequately treated by establishing certain 'moving' jump conditions across 
them. We show how a first integral of the differential equations, if available, can also be used for 
checking the accuracy of the numerical solution. 
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1 Introduction 

When solving a physical problem, one usually encounters a set of coupled nonlinear differential 
equations, called dynamical equations or equations of motion, describing the dynamics of the system. 
This set of differential equations usually emanate from a general physical principle, and might include 
some subsidiary equations which can be categorized in two distinct classes, constraints or integrals 
of motion. The degrees of these subsidiary equations are usually at least one order lower than the 
dynamical equations. At the practical level these two classes of subsidiary conditions are treated very 
differently. The first class has to be solved simultaneously with the rest of the dynamical equations. 
However, the integrals of motion need only be used to put some constraints on the initial conditions. 
Moreover they can be used as a consistency check on the solutions. 

For physically interesting cases, usually time evolution problems appear as initial valued ones, 
and static problems as boundary valued ones. While solving these problems, one often encounters 
various kinds of singularities. These singularities are usually indications of some profound physical 
laws or processes with significant implications. The most common type of singularities are those in 



'Email address: ss-gousheh@cc.sbu.ac.ir 



1 



which the coefficient functions of the differential equations have singularities at some fixed points, 
for example (5-functions. These can be called the fixed type of singularities. In a previous work [1] 
we discussed several general methods to make an efficient numerical algorithm for boundary valued 
problems of this type. 

In this paper we present a method for handling a class of initial valued, coupled, non-linear 
differential equations, whose solutions contain moving singularities. These singularities have the 
property that their positions and severities are apriori unknown and depend on the solutions yet to 
be obtained. Singularities of this type have been noted in such diverse fields as celestial mechanics, in 
particular the classical Kepler problem [2] , and in the study of tensor fields defined on moving surfaces 
[3]. We show that this type of singularity can be adequately treated numerically by establishing 
certain moving jump conditions across them. We establish the accuracy of our numerical solutions 
by showing that the equation representing the integral of motion is satisfied at all values of the 
independent variable including at the positions of the moving singularities. 

For the integration algorithm wc use the basic fourth-order Runge-Kutta method. It is worth 
mentioning that more accurate integration algorithms exist. For example, there are exponential and 
Bessel fitted variable step method of order 6 due to Raptis and Cash [4] , and also a variable step P- 
stable method of order 6 and 8 with a phase lag of the same order due to Simos [5] . For the embedded 
Runge-Kutta, formulae of order 6(5) and 8(7) have been developed by Prince and Dormand [6] and 
formulae of order 8(6) and 12(10) have been developed by Dormand et. al. [7]. These integration 
algorithms should be more efficient for higher accuracy. Our choice of the integration algorithm is 
based on the following reasons. First, our main objective has been to find a solution to the problem 
of moving singularities and not the efficiency of the integration algorithm itself. Second, it turns out 
that in this problem the cumulative error of even the fourth order Runge-Kutta integration algorithm 
is negligible compared to the error introduced to the solutions at each jump across the singularities. 

The set of equations that we discuss results from a classical model of gravitation in Robertson- 
Walker cosmology in which the signature of the metric undergoes a transition from a Euclidean to a 
Lorentzian domain. In section 2 we briefly discuss the physical origin of the problem and show its 
reduction to a set of ordinary differential equations. There, one sees an example in which this set 
automatically includes a subsidiary equation which is an integral of motion, along with the dynamical 
equations. It is worth mentioning that if an integral of motion is not directly included in the set of 
equations of motion, it can some times be derived directly from the dynamical equations. Also we 
employ a reparameterization transformation which allows one to seek continuous solutions across the 
hypersurface of signature change. Moreover, we employ a set of transformations which reduces the 
degree of severity of the moving singularities. The reader who is interested only in the numerical 
methods can skip to section 3, without loss of continuity. 

2 Derivation of dynamical equations 

Traditionally, one of the features of classical gravity is that the signature of the metric is usually 
considered as fixed. If one relaxes this condition, one may find solutions to the field equations which 
exhibit a signature transition [8, 9]. In the model that we study here a real scalar field is taken as the 
matter source interacting with gravity and itself in a Robertson- Walker geometry whose signature 
evolution is controlled by a preferred coordinate. In this model, wc seek solutions to the dynamical 
equations which are smooth and continuous across the hypersurface of signature transition, where the 
metric is degenerate. The alternative would have been to find solutions by solving Einstein's equations 
in disjoint regions next to the hypersurface, and then finding jump conditions to match them [10]. For 



2 



the spatially flat universes, the first approach yields exactly solvable Einstein's equations [8]. Here, 
we discuss the general case which includes the spatially flat as well as non-flat cases and solve the 
resulting dynamical equations numerically. For more details of the physical basis and signiflcance of 
the problem, we refer the interested reader to the reference [9]. 

Consider gravity coupling to a scalar field through Einstein's equation, 

Giiv = K^^i/M, (1) 
where the scalar field ^ is a solution of the Klein-Gordon equation, 

dU 

A0-^=O. (2) 

Here, 0^^,^ is the Einstein tensor constructed from torsion-free connections compatible with the metric, 
and U{(j)) is the scalar potential for the real scalar field ^, which interacts with itself and gravity 
through the stress-energy tensor T[(j)\. 

The above coupled equations are to be solved in a domain that would lead to Robertson- Walker 
cosmologies with Lorentzian signature. However, if the metric is suitably parametrized, one expects 
to see continuous transition to a Euclidean domain. As in [8], we adopt a chart with coordinate 
functions } where the hypersurface of signature change would be located at /? = 0. The 

metric can be parametrized to take the form 

g = -Pdp ^d(3+ J^jj^^^y^, E dx' dx\ (3) 

where = We seek solutions of the form R = R{f5) and (j) = (l){(3). Now, it is apparent that 

the sign of f3 determines the geometry, being Lorentzian if /? > and Euclidean if /? < 0. For f3 > 0, 
the traditional cosmic time can be recovered by the substitution t = (2/3)/3^/^. Adopting the chart 
{t, x'} and using equations (1) through (3) with units in which k = 1, one finds 

-3|,-3^ + \ + U{ct.)=0, (4) 

* + 4 + f = W 

where a dot represents differentiation with respect to t and — oo<0<oo,O<i?<oo. Now a 
solution to the problem is furnished by finding R{t) and (p{t), for a given U{(p). Note that these 
equations are not all independent. For example equation (5) can be obtained by combining equations 
(4) and (6). Upon a closer inspection we recognize that this is due to the fact that equation (4) 
is not a dynamical equation, rather it is actually an integral of motion representing a zero energy 
condition. That is, any solution of the dynamical equations (5) and (6) would yield a constant total 
energy (equation (4)). However, Einstein's equations demand zero energy solutions only. 

As is apparent from the dynamical equations, we have moving singularities at all times for which 
R = 0. These moving singularities are potentially very severe and, as we shall see later, (/) actually 
diverges there. Wc can get an indication on the divergence of (p from equations (4,5). These equations 
indicate that U{(p) has to cancel the divergence of the k/R^ terms, and for all physically relevant 
potentials, this implies that cf) has to diverge. We therefore need to use a set of transformations to 
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reduce the severity of the divergence of solutions. We expect the following transformations to render 
the solutions more manageable, since it is formed of products of factors which go to zero and infinity 
at about the same strength, 

X = i2^/^cosh(a0), (7) 
Y = R^/'^smh{a(j)), (8) 

where = |. 

The above equations are considerably simplified if we take the potential to be 

2a^{X^ - Y^)U{cl){X, Y)) = aiX^ + asF^ + 2bXY, (9) 

where 04 , 0,2 and b are adjustable parameters. This choice for the potential stems from the fact 
that the left hand side of equation (9) directly appears in the Lagrangian, from which the dynamical 
eqations can be derived. The features of this potential and the physics involved in the choice of its 
parameters have been discussed in [8, 9]. 

The dynamical equations (5) and (6) in terms of X and Y and the evolution variable /3 now 
become, 

Y" = \ (^) Y' - \pkY{X^ - y2)-2/3 _ p^^^Y + (10) 

X" = \ X' - ^(3kX{X' - Y^r^l^ + (5{a,X + bY), (11) 

subject to the subsidary 'zero energy condition', equation (4), which can be written in terms of the 
new variables as 

i-X'' + Y'') - ^k{X'^ - y2)V3 + („^j^2 ^ ^^y2 ^ 2bXY) = 0. (12) 

Here, a prime represents differentiation with respect to (3. Note that this equation does not contain 
any singularity, and equations (10) and (11) are actually less singular than equations (5) and (6). 
The coupled equations (10) and (11) must now be solved and, as explained before, equation (12) is 
merely a restriction on the initial conditions. However, it can also be used as a consistency check on 
the analytical or numerical solutions. These equations do not seem to have a closed form solution so 
a numerical treatment is necessary. 

3 The numerical method 

The dynamical equations that we have to solve are equations (10, 11). As mentioned in the last 
section, equation (12) is an integral of motion. That is, any true solution to equations (10, 11) 
automatically satisfies equation (12) for all values of the independent parameter, if it satisfies it at 
any one point. Therefore, we can use equation (12) to put a restriction on the initial conditions. 
More importantly, one can check the accuracy of the solutions by seeing how well equation (12) is 
satisfied, as the algorithm integrates the dynamical equations. 

As a first step towards a numerical solution to the equations, we should study the restrictions 
imposed by the set of differential equations on the initial conditions. These restrictions are the result 
of the requirement of consistency of the initial conditions with the dynamical equations. However we 
can accomplish a more complete task by finding the general form of the analytic solutions close to 
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the initial point. These solutions certainly include the complete information on the allowed set of the 
initial conditions^. Moreover, the knowledge on the analytic solutions help with the first few steps of 
the integration algorithm. 



3.1 Analytic solutions close to the initial point 



In order to find analytic solutions which are valid near the initial point (/? = 0), we first study the 
restrictions imposed by equations (10)-(12) on the initial conditions. This is done by noting that in 
order to have well behaved solutions close to /? = 0, the first term of equation (12) shows that we must 
either have ~ and Y'{(3) ~ where n^, Uy > 1/2, or \X'{0)\ = \Y'{0)\. However, the 

first terms on the right hand sides of equations (10) and (11) impose a more severe restriction. These 
two equations admit solutions X'{P) ~ P^/^ and Y'{/3) ~ /3^/^ close to /3 = 0, however, this class of 
solutions does not admit real or solutions across (3 = 0. One can show that regular solutions close 
to /3 = are of the form 



XiP) = A^p' + Xo, where = - 



kXn 



AyP^ + Yo, where Ay 
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4 {Xi - 1-0^)2/3 

3 kYp 
"4(X2-ro')2/3 



+ aiXo + 6^0 
- a2lo - ^-'^o 



(13) 
(14) 



with Xq = X{0), etc. Therefore, the initial conditions on the first and second derivatives must satisfy 
the relations 



X'{0) = y'(0) = and X"{0) = Y"{0) = 0. 



(15) 



Strictly speaking the conditions on the second derivatives are not initial conditions but rather 
consistency checks, since we have coupled second order equations. Therefore, the initial values for 
the functions X and Y must now satisfy, c.f. equation (12), 



Y^)^f^ + (aiX2 + as^o' + 26X0^0) = 0. 



(16) 



The contour plots of equation (16) for k = ±1 arc given in figure 1. Along the contours, one finds the 
possible initial values for X and Y. Although equation (16) is equivalent to a sixth order algebraic 
equation which cannot be directly solved analytically, we can solve it by going back to the original 
variables R and 0. The solutions are either R{0) = giving (/)(0) = ±00, which we exclude because 
we have been seeking continuous solutions across /3 = 0, or R{0) / (it is a free parameter) with 



^(0) = ^ cosh"^ 
2a 



(17) 



where 



D 



9k 



ai — a2 



4R{0y 



and B 



m 
T 



Therefore, the acceptable values of Xq and Yq ( -'^o > |^o|) can also be obtained analytically from 
equation (17). 

'^The question of the allowed set of the initial conditions, though interesting enough in its own right in all problems 
of this type, is of crucial importance for the problem at hand, as the determination of the correct initial conditions is 
an open problem in cosmology. 
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Figure 1: The contour plots of the allowed initial values of X and Y , satisfying the equation of constraint (16) for 
k = ±1. The point (0,0) is a solution and the curves approaching this point actually pass through it, although this is 
not shown on the plots due to the limitations on the numerical accuracy. 



3.2 Integration algorithm 

The important feature of equations (10) and (11) is that they are singular for all (3 at which X = . 
At these critical values of (3 (Pc), the original variables take the values R{l3c) = and 0(/3c) = ±oo. 
We can directly infer from the differential equations that the solutions for the new variables and 
their first derivatives have to be continuous across the singularities. However, the second and higher 
derivatives will be singular at Z?,.. That is, the singularities of the new variables are considerably 
milder than those of the original variables. 

Although the solutions and their first derivatives are continuous across Pc, they cause problems 
for the integration algorithm. Any attempt in solving these equations involves handling these moving 
singularities, as one encounters them when integrating the coupled equations. To proceed, wc first 
establish jump conditions across these singular points as follows: close to (3c we assume that the 
solutions have the following linear forms 

X± = a± + b±l3, (18) 
Y± = c± + d±p, (19) 

where it refers to the right or left hand sides of the singularity, respectively. Substituting the above 
equations in (10) and (11) and dropping all non-singular terms, one can integrate these equations 
in the interval Pc ~ Pc + where 2e is the distance across the jump. For the integration we have 
dropped all terms which would give rise to contributions 0{e^^^). One finds at Yc = ±Xc 



9 (2Xe6)^/^ 
■4 (6_^d_)2/3' 



b+-b^ = --k-^^-^—rPc = Tid+-d-), (20) 



where e can be taken as small a value as is desired for any required accuracy. Equation (20), together 
with the requirement of continuity of X and Y, establish our jump condition for handling the sin- 
gularities of the differential equations. It is apparent from equation (20) that the slopes X'{P) and 
Y'(P) are continuous at /Re- 
writing an actual algorithm for handling these singularities requires some care. Let us first write 
the original variable R{P) in terms of the new variables. 



R = {x^ - y^) 



2x1/3 
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Figure 2: Solutions for X{/3) (broken curve) and Y{/3) (solid curve), for k = ±1. The values of the parameters are 
6 = 2, A = 0, = 4.5. 



We recall that the differential equations become singular when R = 0. As we integrate them, when 
\R\ becomes small (less than 1) we reduce the step size by one order of magnitude since the crossing of 
R{f3) through zero at (3c is rather steep. Then at the first instant when \R\ becomes smaller than 0.1, 
henceforth called 'the fixed point', the algorithm records all the relevant values {X,Y, X' ,Y' , R, (3) 
and continues integrating towards the singular point with yet finer steps. Past the fixed point, the 
singular terms in the differential equations become too large and no integration algorithm can give 
reliable values for X and Y. However, we can use the information obtained past this point to pinpoint 
Pc as follows: We record the last two values of (3 and R right before the instant when the sign of 
R changes. Then assuming R (x {P — PcY^'^ (which is consistent with equations (18) and (19)) one 
obtains 

Rlf3i - Rll32 

Rl-Rl ■ ^ > 

Having obtained /3c, we can calculate e = /3fixed ~ Pc and use a linear extrapolation to obtain Xc and 
Yc and see whether Yc = =bXc as a consistency check. We can then calculate the values of the slopes 
{X\ Y') on the other side of the singularity (at /3over = /^c + e) using the jump conditions (equation 
(20)) and then using linear extrapolation on both sides of the singularity, the values of the functions 
can be calculated at /3over- The algorithm then continues integrating with fine steps until the value 
of increases beyond 1 and with regular steps until it approaches the next singularity. 

We use a set of parameters (6 = 2, A = 0, m? = 4.5) in equations (10-12) which arc physically 
relevant and choose our initial conditions consistent with equations (15) and (16). Recall that since 
equation (12) is a constant of motion, if it is satisfied at /? = 0, for a true solution it will be satisfied 
at all other values of (3. Therefore if the initial values for the functions X and Y satisfy equation (16) 
at /3 = 0, equation (12) should always be satisfied. For integrating equations (10) and (11), we have 
used the fourth order Runge-Kutta method. The resulting solutions for k = ±1 arc shown in figure 

2. It is apparent from figure 2 that at the singular points X = zty, the solutions are continuous and 
the singularities are very mild. As a measure of the accuracy of the solutions, we have computed 
the 'zero energy condition', equation (12), as a function of (3 for k = ±1 which are shown in figure 

3. It is evident from figure 3 that the values of 'total energy' stay very close to zero, thus indicating 
the validity of the numerical solution. In figure 4 the variations of the original variables (f) and R 
are shown. As can be seen, ^ actually diverges at the singular points. As a further check, we have 
numerically recovered the analytic solutions presented in [8] for A; = in every detail. 
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Figure 3: The graph of the total energy defined by equation (12) for k = ±1. As is apparent from the graphs, the 
zero energy condition is satisfied to a high accuracy. The small jumps in the graphs are at the critical values of (3 where 
there are singular points: Y{f3c) = ±X(/3c)- 
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Figure 4: Graphical representation of the original variables (piP) and R{0), for 6 = 2, A = 0, = 4.5, and k = 1 
(broken curves) and fc = — 1 (solid curves). 
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4 Conclusions 



We have shown how a particular class of initial valued coupled non-linear ordinary differential equa- 
tions with moving singularities can be numerically solved. The main obstacle of having moving 
singularities can be overcome by establishing a set of jump conditions across them. These conditions 
are obtained by approximating the form of the solutions close to the singular points and directly 
integrating the differential equations in the neighbourhood of these points. We have found that a 
first order approximation close to these points is sufficiently accurate. Also since the main source 
of error in the solutions eminates from the jump conditions, we have found that the fourth order 
Runge-Kutta is sufficient for the integration algorithm. 
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